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The methods of density-functional perturbation theory may be used to calculate various physical 
response properties of insulating crystals including elastic, dielectric, Born charge, and piezoelectric 
tensors. These and other important tensors may be defined as second derivatives of the total energy 
with respect to atomic-displacement, electric- field, or strain perturbations, or as mixed derivatives 
with respect to two of these perturbations. The resulting tensor quantities tend to be coupled in 
complex ways in polar crystals, giving rise to a variety of variant definitions. For example, it is gen- 
erally necessary to distinguish between elastic tensors defined under different electrostatic boundary 
conditions, and between dielectric tensors defined under different elastic boundary conditions. Here, 
we describe an approach for computing all of these various response tensors in a unified and system- 
atic fashion. Applications are presented for two materials, wurtzite ZnO and rhombohedral BaTi03, 
at zero temperature. 
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I. INTRODUCTION 

The methods of density- functional theory (DFT) 1 and 
density-functional perturbation theory (DFPT) 2 ' 3 have 
been shown to give a successful description of the dielec- 
tric and piezoelectric properties of a wide range of mate- 
rials in which electronic correlations are not too strong. ,5 
Many properties of interest can be computed directly 
from DFT using finite differences - for example, elas- 
tic constants computed from the stress arising from a 
small applied strain, or dynamical effective charges com- 
puted from polarizations 6 arising from small sublattice 
displacements. On the other hand, the use of DFPT 
methods is becoming increasingly popular because it can 
be used to compute such response properties directly, 
without the need for multiple ground-state calculations, 
thus providing the desired response properties in a more 
automated, systematic, and reliable fashion. 

As a result, improved DFPT capabilities have been 
implemented in recent years in several of the com- 
puter code packages commonly used by the computa- 
tional electronic-structure community. 7-9 This develop- 
ment has been most thorough in the case of the open- 
source ABINIT computer package, 7 in which the capa- 
bility for handling strain perturbations 10 has recently 
been added to the previous implementation of atomic- 
displacement and electric-field perturbations. This de- 
velopment opens the prospect for a systematic treatment 
of three kinds of perturbations in insulating crystals on 
an equal footing: periodicity-preserving atomic displace- 
ments (i.e., zone-center phonons), homogeneous electric 
fields, and homogeneous strains. These three degrees of 
freedom are often strongly coupled, especially in polar 
materials used in modern ferroelectric, piezoelectric, and 
dielectric applications. 



In this work, we show that such a systematic approach 
is now not only practical, but especially powerful. That 
is, we show that computing the full set of six elemen- 
tary (or "bare") response tensors (force-constant, dielec- 
tric, elastic-constant, Born-charge, internal-strain, and 
piezoelectric tensors) associated with these three kinds 
of perturbations provides an extremely valuable database 
that can subsequently be used to compute a wide variety 
of relevant physical properties. Among these, for exam- 
ple, are the physical dielectric, elastic, and piezoelectric 
tensors (in which atomic displacements have been taken 
into account), elastic compliances, free-stress dielectric 
tensors, fixed-electric-displacement elastic tensors, alter- 
native piezoelectric tensors, and electromechanical cou- 
pling constants. The ability to access this wide range of 
secondary properties becomes possible only after the full 
set elementary response tensors has been systematically 
computed. 

The rest of this paper is organized as follows. In Sec. II, 
we present the formal development, denning the vari- 
ous elementary response tensors and showing how other 
response tensors of interest can be derived from these. 
Then in Sees. III-IV we apply our approach to wurtzite 
ZnO and rhombohedral BaTiOa as two paradigmatic ex- 
ample systems. We first briefly describe the practical 
details of the calculations in Sec. Ill, and then present 
the results for the ground-state properties, elementary 
response properties, and derived response properties, in 
Sec. IV. We conclude with a summary in Sec. V. A care- 
ful formulation of the theory for the case in which strains 
and electric fields are simultaneously present is deferred 
to the Appendix. 
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II. FORMALISM 
A. Elementary response tensors 

Consider an insulating crystal with N atoms per unit 
cell. We choose a reference state in which the lattice 
vectors are ai, &2, and a3, the cell volume is Qq, and the 
atomic coordinates are Rm^ . Here m is a composite label 
(atom and displacement direction) running over 1, 37V, 
and we assume that this structure is the equilibrium one 
in vanishing macroscopic electric field. 

We consider three kinds of perturbations applied to 
such a crystal: (i) displacements u m of the atoms away 
from their equilibrium positions, (ii) homogeneous strains 
rjj where j = {1...6} in Voigt notation, and (iii) homoge- 
neous electric fields £ a where a = {x, y, z} are Cartesian 
directions. We restrict our discussion to atomic displace- 
ments that preserve the primitive-cell periodicity, i.e., to 
degrees of freedom corresponding to zone-center phonon 
modes only. Also, we will restrict ourselves entirely to 
zero-tempcrature properties. 

The corresponding responses that are conjugate to 
these three perturbations are (i) forces F m , (ii) stresses 
<jj, and (iii) polarizations P a . From these, one can con- 
struct the response functions of primary interest: "diag- 
onal" responses K mn = dF m /du n (force-constant ma- 
trix), Xap = dP a /d£p (dielectric susceptibility), and 
Cjk = d<jj/dr]k (elastic constants), and "off-diagonal" re- 
sponse tensors Z ma = dP a /du m (Born effective charge), 
A m j = dF m /drjj (internal strain), and e a j — dP a /drjj 
(piezoelectric response). However, in order to define 
these quantities carefully, it is important to clarify the 
constraints or boundary conditions that apply to each 
definition. For example, the elastic constants Cjk may 
be defined allowing or not allowing internal atomic dis- 
placements ( "relaxed- ion" or "frozen- ion" ) , or under con- 
ditions of fixed electric (£) or displacement (D) field. 

We take the approach here of systematically defining 
all response properties as appropriate second derivatives 
of the energy E per unit volume with respect to the 
perturbations. To be more precise, in the presence of 
strains we define E as the energy per undeformed unit 
cell volume fio, while in the presence of electric fields E 
is modified to become an electric enthalpy 11 by adding 
a term proportional to — P • £, where P is the electric 
polarization. 6 (While a direct treatment of finite ^-fields 
is now possible, 12,13 only infinitesimal ^-fields need to be 
considered here.) In general, we define E as 



E(u,£,r]) = -|- 



(1) 



where i^eii * s ^ ne usua l zero-field Kohn-Sham energy per 
cell 14 of the occupied Bloch functions and ft is the de- 
formed cell volume. However, when strains and electric 
fields are simultaneously present, care is needed in the 
interpretation of Eq. (1); this is explained in the Ap- 
pendix, where a more precise formulation is given in the 



form of Eq. (All), which supersedes Eq. (1). In short, 
the difficulty is connected with the distinction between 
"proper" and "improper" piezoelectric constants; 15 we 
should like our formulation to lead to the former and not 
the latter. The factor of O/rio has been inserted in the 
last term of Eq. (1) towards this purpose, but this is not 
sufficient by itself. In addition, Eq. (1) should be rewrit- 
ten in terms of "natural variables" w, and r/, where 
£ has been replaced by a reduced electric field £' that is 
defined in Eq. (A5). When partial derivatives are taken 
with respect to these natural variables, one automati- 
cally obtains the "proper" piezoelectric tensors. Indeed, 
as explained in the Appendix, all appearances of £ should 
be replaced by £', with a similar replacement for polar- 
izations, in the remainder of this paper. However, for 
the sake of clarity of presentation, this notation has been 
suppressed in the main body of the paper. 

Accordingly, we provisionally write E = E(u 7 £,r/) as 
a function of arguments u m , £ a , and rjj, with the under- 
standing that the notation of the Appendix supersedes 
the notation used here whenever strains and electric fields 
are simultaneously present. We then expand around a 
zero-field reference system as 



E : 



+ \B. 
+ B. 



-f- A a £ a + Aj T}j 



\Bj k rjj ri k 



ma. ""in <--a 



+ terms of third and higher order . 



(2) 



We use an implied-sum notation throughout. In this ex- 
pansion, the first-order coefficients A m , A a , and A } en- 
code the forces (F m = — floA m ), polarizations (P a — 
— A a ), and stresses (aj — Aj), respectively. (Hence- 
forth we shall assume that the atomic coordinates and 
strains are fully relaxed in the reference system, so that 
A m = Aj = 0.) The diagonal- block second-order coeffi- 
cients B mni B a /3 and Bjk and off-diagonal second-order 
coefficients B ma , B m j , and B a j correspond to the force- 
constant, elastic-constant, and susceptibility tensors, and 
to the Born-charge, internal-displacement, and piezoelec- 
tric tensors, respectively. 

Inserting appropriate signs and cell- volume factors, the 
elementary second-derivative response-function tensors 
arc defined as follows. The force-constant matrix 



K mn — fio 



d 2 E 



du m du n 

the frozen-ion dielectric susceptibility 

8 2 E 



Xaf3 = - 



d£ a d£p 

and the frozen-ion elastic tensor 



U,7] 



drjjdrik 



(3) 



(4) 



(5) 
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are the elementary diagonal-block tensors, while the off- 
diagonal blocks are the Born dynamical effective charge 
tensor 



Z m a — — ^0 



d 2 E 

du rn d£ a 



(6) 



the force-response internal-strain tensor 



A 



-IV 



d 2 E 



du m dr]j 

and the frozen-ion piezoelectric tensor 



(7) 



d 2 E 

dE a drjj 



(8) 



The bar on quantities Xafi-, Cjk, and e a j indicates a 
frozen-ion quantity i.e., the fact that atomic coordinates 
are not allowed to relax as the electric field or homoge- 
neous strain is applied. Note that the clamped-ion elas- 
tic tensor Cjk and piezoelectric tensor e a j are generally 
not physically relevant quantities, except in cases of high 
symmetry where atomic displacements do not occur to 
first order in strain. The clamped-ion susceptibility ten- 
sor Xap is the purely electronic one that is measured in 
response to AC or optical fields at frequencies well above 
the phonon-frequency range (corresponding to e°° in the 
polariton language). 

The force-response internal-strain tensor A m j must be 
distinguished from the displacement-response internal- 
strain tensor T n j = A m j(i4T _1 ) mn that describes the first- 
order displacements resulting from a first-order strain; 
both occur in the literature, frequently without careful 
differentiation. The piezoelectric tensor e a j (often de- 
noted alternatively as c a j) describes the change of po- 
larization arising from a strain, or a stress arising from 
a change of £-field, while the d, g, and h piezoelectric 
tensors are defined under different constraints and have 
slightly different physical meanings. 16 Finally, we remind 
the reader that there is some subtlety in the definition 
of the piezoelectric tensors related to the specification of 
the energy functional when both fields and strains are 
present, leading to a distinction between "proper" and 
"improper" piezoelectric constants 15 as will be discussed 
more fully in the Appendix. Throughout this paper, we 
adopt the convention that all piezoelectric tensors are 
"proper" ones unless otherwise noted. 

We shall refer to the quantities defined in Eqs. (3-8) 
as the "elementary" or "bare" response tensors. These 
are the quantities that will be calculated once and for all 
using the DFPT capabilities of a code package such as 
ABINIT. All of the derived tensor properties described 
in the following subsections can then be calculated from 
these using simple matrix manipulations, as we shall see. 



B. Relaxed-ion tensors 

Generally, the physical static response properties of 
interest must take into account the relaxations of the 
ionic coordinates. This becomes especially important for 
non-centrosymmetric polar systems, such as ferroelectric 
ones, where these various effects become coupled. In- 
stead of "clamped-ion" quantities x, C and e defined at 
fixed u, we can define the corresponding "relaxed-ion" or 
"dressed" response tensors C, x 5 an d e as follows. To 
develop expressions for these, we let 



E(t],£) = min E(u,£,r)) . 



(9) 



Referring back to Eq. (2), setting dE/du n = 0, 
dE/d£ a = 0, and dE/dr]j = 0, and assuming that the 
reference configuration is one in which the forces A m van- 
ish, we find 

B nrn u m + B na £ a + E n j Tjj 
from which it follows that 

u m = —{B~) mn [B nj tjj + B na £ a ] . (10) 



Defining 



Xa/3 



d 2 E 

d£ a d£fj 

d 2 E 
df]jdr] k 

8 2 E 



dS a dr}j 



(ii) 

(12) 
(13) 



and using Eqs. (3-8), we find that the physical relaxed- 
ion dielectric, elastic, and piezoelectric tensors become 



Xck/3 — "Ka.fi H~ ^ma. )mn ^n(3 •> 

Cjk ~ A m j (K )mn A-nk •> 



(14) 
(15) 
(16) 



respectively. 

Note that Eqs. (14-16) cannot be naively evaluated as 
written because the force-constant matrix K is singu- 
lar, due to the fact that K has three vanishing eigenval- 
ues associated with translational symmetry. Moreover, 
in soft-mode systems, other eigenvalues may be close to 
zero, and care should be taken to distinguish these from 
the translational ones. For these reasons, we have imple- 
mented a careful procedure for obtaining the "pseudo- 
inverse" of K; throughout these notes, whenever we refer 
to K , we really mean the pseudo-inverse. 

We proceed as follows, (i) We identify the three- 
dimensional space of acoustic modes (i.e., uniform trans- 
lations), and construct a (3N) x (3JV) orthogonal matrix 
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U whose first three columns correspond to these trans- 
lational modes; the remaining columns of U are formed 
by Graham-Schmidt orthogonalization of the basis, (ii) 
We construct K 1 = U KU T , whose upper 3x3 block rep- 
resents the acoustic subspace and therefore ought to be 
zero, (iii) We let K' rcd be the lower (37V - 3) x (3N - 3) 
block of K' , corresponding to the reduction to the com- 
plementary subspace of optical modes, (iv) We invert 



Following a line of reasoning similar to that leading from 
Eq. (9) to Eqs. (14-16) and setting <jj — 0, one obtains 



(°0 

X a f3 



Xafi + e a j (C 1 ) J k epk ■ 



(20) 



Typically, an AC dielectric measurement will access the 
true static susceptibility x^ as long as the frequency is 
much less than that of sample resonances (elongational, 



-^-rcd by standard means (taking appropriate measures in bending, or torsional modes), and x^ a t frequencies 

much higher than sample resonances (but much less than 
phonon frequencies). 

Before leaving this subsection, we note that it is con- 
venient to define inverse dielectric tensors 



case this matrix is nearly singular, as when soft modes 
have nearly vanishing frequencies). Let the result be de- 
noted as (-K" _1 )£ ed . (v) We pad (if -1 )^ with zeros in 
the first three rows and columns to form the (3iV) x (37V) 
matrix (TV -1 )', (vi) Finally, we define the pseudo-inverse 
of K to he K- 1 =U T (K- 1 )'U. 

Thus, by construction, the resulting pseudo-inverse 
K^ 1 is zero in the subspace of translational modes, and 
is the inverse of the original matrix in the complemen- 
tary subspace. As a result, any time K^ 1 is multiplied 
by another tensor, a pre-projection onto the complemen- 
tary subspace of dimension 3iV — 3 is effectively carried 
out. In other words, the acoustic sum rule is effectively 
enforced in any operation involving K~ x . 



C. Other derived tensor quantities 

In the previous subsection, we showed how to obtain 
the static dielectric susceptibility tensor Xa/3, the elastic 
tensor Cjk, and the piezoelectric tensor e a j. These quan- 
tities are defined under conditions of controlled strain 
and electric field. From these three ingredients, it is 
straightforward to form many other useful tensor quan- 
tities describing physical properties defined under other 
constraints or boundary conditions, as we shall see in this 
section. 



1. Dielectric tensors 

The susceptibility tensor Xap is defined at fixed (van- 
ishing) strain; the corresponding dielectric tensor is 



£0 (Sa/3 + Xap) 



(17) 



where e is the susceptibility of free space (SI units are 
used throughout) and the superscript (77) indicates that 
the derivative is taken at fixed strain. Often one is inter- 
ested instead in the free-stress dielectric tensor 



(18) 



which incorporates the free-stress susceptibility x^ ■ An 
expression for the latter is easily derived from the elastic 
enthalpy 

H(a,£)=min[E(r ] ,£)-r ]j a j ] (19) 



(e W)-i 



(21) 
(22) 



for later use. 



2. Elastic and compliance tensors 

The elastic tensor Cjk defined in Sec. II B is the one de- 
fined under conditions of fixed (vanishing) electric field: 

Cjk = C-J. It may sometimes be of interest to treat 



instead the elastic tensor C-j^ defined under conditions 
of fixed electric displacement field D. For example, in 
the case of a thin film of dielectric material sandwiched 
between much thicker layers of other insulating host ma- 
terials, the electrostatic boundary conditions fix the com- 
ponent of D, not £, normal to the interfaces. One readily 
obtains 



(D) 



c jk + e "i P a p e /?fe 



(v) 



(23) 



It is also straightforward to obtain the corresponding 
elastic compliance tensors either under zero f -field 



or under zero D-ficld 
boundary conditions. 



(24) 



(25) 



3. Piezoelectric tensors 

The formulation of an energy functional appropriate to 
the simultaneous treatment of strains and electric fields 
is rather subtle, as discussed in the Appendix. There, we 
show that the proper (relaxed-ion) piezoelectric tensor 
Cjk introduced in Sec. II B may be written as 



dP n 



dr)j 



(26) 
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or equivalently, by a thermodynamic relation, 17 ' 16 



d£ n 



(27) 



where it is understood (see Appendix) that in these and 
subsequent equations, P and £ are to be interpreted as 
the reduced polarization P' and electric field £' and of 
Eqs. (A4) and (A5), respectively. This is done so that 
Eqs. (26) and (27) will yield the proper, rather than the 
improper, piezoelectric tensor. 15 

In view of Eq. (27), e a j is sometimes referred to as 
a "piezoelectric stress constant." In any case, it is the 
natural piezoelectric constant defined under conditions 
of controlled £ and rj. On the other hand, the "piezo- 
electric strain constant" d a j , defined under conditions of 
controlled £ and er, is equally or even more commonly 
discussed in the literature; it is defined via 



d, 



or 



drjj_ 
d£ a 

dP a 



d<jj 



and is given in terms of e via 

A - S {£} p , 

u aj — jk a & 



(28) 



(29) 



(30) 



as can again be derived from thermodynamic rela- 
tions. 17,16 Two other, less commonly used, piezoelectric 
tensors are g a j and h a j , defined under conditions of fixed 
(D,a) and (D,r]), respectively, and given by 



9aj 



Pa0 d P3 



(31) 
(32) 



where the (3 are the inverse dielectric tensors defined in 
Eqs. (21-22). These have the properties 

fyj = 9aj 5D a , 
S£ a = —g a j 5a j , 



-h a j 5D C 



5£ a = —h a j Sr/j 



(33) 



label j. For example, is a dimensionless measure 
of the coupling of electric and strain degrees of freedom 
along the z axis. Roughly speaking, a coupling factor 
close to unity implies an excellent impedance match for 
the material used as an electromechanical transducer be- 
tween the specified electric and elastic channels (a cou- 
pling factor greater than unity is forbidden by stability 
considerations) . 16 

Note that k a j in Eq. (34) does not transform like a 
tensor, and the usual implied sum notation does not ap- 
ply to this equation. Instead, we can define a tensorially 
correct, dimensionless coupling tensor via 



JC = [/?W]V2 .d.[C< £ >]V2 , 



(35) 



where an obvious matrix-product notation is used. The 
standard "singular-value decomposition" can be used to 
write K as 



<k x 0^ 

a: u I k , o 

k 3 0/ 



■V 1 



(36) 



defined by requiring that U and V be orthogonal 3x3 
and 6x6 matrices respectively, and that the k u should 
be positive. (Alternatively, the k~ can be determined 
as the eigenvalues of the 3x3 symmetric matrix K, IC T = 
[pW] 1 / 2 dC& (f [Z^)] 1 / 2 .) For each singular value, the 
corresponding columns of U and V give the pattern of 
electric field and of strain, respectively, that are directly 
coupled to one another by JC. 

The coupling factors can be related to differences be- 
tween dielectric or compliance tensors defined under dif- 
ferent boundary conditions. Starting from Eqs. (18-20) 
and Eqs. (23-25), one can show 

e («x) _ e (v) = dC (£) d T = [ c M]i/2 jcjcT [e W]i/2 ; (37) 



S (£) _ S (D) = d Tp(.) d = [s ( £)] l/2 K T K 



/2 



(38) 



Specializing to high-symmetry situations in which e is 
necessarily diagonal, one finds, for example, 



>) 

tota 



taa 



M 

taa _ 72 



(39) 



4- Piezoelectric coupling coefficients 

The most common definition of the piezoelectric cou- 
pling factor k a jis given by 16,18 




This applies to the case where the field is applied only 
along a and the only non-zero stress is the one with Voigt 



III. METHODS AND DETAILS OF THE 
CALCULATIONS 

Our ab-initio calculations were carried out using the 
ABINIT code package. 7 ' 19 Specifically, we first carried 
out full structural relaxations for both materials. Next, 
response-function calculations were carried out in order 
to obtain first derivatives of the occupied wavefunctions 
with respect to the perturbations of atomic displace- 
ments (i.e., phonons at q = 0), uniform electric field, and 
strain. These were then used to compute the elementary 
second-derivative response-function tensors, Eqs. (3-8), 
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of Sec. II A. Except for the diagonal elements of some el- 
ementary tensors, this was done using a non-variational 
expression that only requires input of one of the two cor- 
responding wavefunction derivatives. 2 (As for the "mixed 
derivative" tensors A, e and Z, strain derivatives were 
used for A and e, and displacement derivatives were used 
for Z. 20 ) Finally, from these elementary response tensors, 
the various secondary response tensors of Sees. II B and 
II C are obtained according to the formulas given there. 
All calculations are at zero temperature. 

The DFT and DFPT calculations for ZnO and 
BaTiOa were carried out using Troullicr-Martins 
pseudopotentials 21 and a plane-wave energy cutoff of 50 
hartree. The Zn pscudopotential includes the 3d elec- 
trons in the valence, as this has been shown to be impor- 
tant for accurate results. 22 An 8x8x8 and 6x6x6 
Brillouin-zone k-point sampling were used for ZnO and 
BaTiC>3 respectively. The exchange and correlation ef- 
fects were treated within the local-density approximation 
(LDA) in the Ceperley- Alder 23 form with the Perdew- 
Wang 24 parameterization. 

Finally, we made one additional modification in the 
case of BaTiOa, where it is well known that the usual 
underestimation of the equilibrium lattice constant asso- 
ciated with the local-density approximation has an un- 
usually profound influence on the ferroelectric distortion, 
which is very sensitive to cell volume. 25 Therefore, to 
get more physically meaningful results, we carried out 
the initial structural relaxation with the cell volume con- 
strained to be that of the experimental structure at zero 
temperature. This is similar in spirit to the use of a 
"negative fictitious pressure" that is a standard feature 
of many first-principles based studies of ferroelectric per- 
ovskite materials. 25 



IV. RESULTS FOR TWO SAMPLE SYSTEMS: 
ZnO AND BaTiOa 

In this section, we consider two paradigmatic systems, 
wurtzite ZnO and rhombohcdral BaTi03. The elec- 
tromechanical properties of ZnO make it a widely used 
material in mechanical actuators and piezoelectric sensor 
applications, while BaTi03 is a prototypical perovskitc 
ferroelectric material. It is of particular interest to com- 
pare and contrast the behavior of these two materials in 
view of the fact that BaTi03 is a soft-mode system, while 
ZnO is not. This may help provide insight into the role of 
the soft mode, which can be expected to lead to enhanced 
piezoelectric and dielectric response and enhanced cou- 
plings. We first describe the results of our ground-state 
DFT calculations, and then present the results for the 
various linear-response tensors as defined in Sec. II. 

Because ZnO is not a soft-mode system, its properties 
depend only weakly on temperature, so that it is not un- 
reasonable to compare room-temperature experimental 
results with zero-temperature theory. BaTi03 is a diffcr- 



TABLE I. Structural parameters of ZnO. Lattice constants 
a and c in A; u is dimensionless internal parameter. 





a 


c 


c/a 


u 


Present work 


3.197 


5.166 


1.616 


0.380 


Previous theory a 


3.286 


5.241 


1.595 


0.383 


Expt. b 


3.250 


5.210 


1.602 


0.382 



a Ref. 28. 
b Ref. 29. 



ent matter, as its properties depend crucially on tempera- 
ture. The room-temperature tetragonal phase has indeed 
been thoroughly studied, and as a result, its dielectric, 
elastic, and piezoelectric constants have been systemati- 
cally measured. 26,27 However, there are formidable diffi- 
culties associated with preparing single-crystal, single- 
domain samples of the low-temperature rhombohcdral 
phase, and of carrying out dielectric and elastic measure- 
ments on such samples at low temperature under well- 
defined electric and elastic boundary conditions. As a 
result, almost no reliable experimental values are avail- 
able for the corresponding materials constants at very low 
temperature. Therefore, for the purposes of this paper, 
we have adopted the approach of providing comparisons 
with experiment for ZnO wherever possible to benchmark 
our approach, and of presenting our calculations of the 
low-temperature properties of BaTi03 as predictions for 
a system that is difficult to characterize experimentally. 

A. Relaxed structural properties 

1. ZnO 

The ground state of ZnO is a tetrahedrally coordinated 
wurtzite structure (space group P63ITIC, point group Cq v ) 
with four atoms per unit cell. The structure is deter- 
mined by three parameters: the hexagonal edge a, the 
height of the prism c, and the internal parameter u. The 
structural results from our full relaxation are given in 
Table I. For comparison, an ideal wurtzite with exactly 
tetrahedral angles an d eq ual-length bonds would have 
u = 3/8 and c/a = ^/8/3. As is typical of DFT calcu- 
lations, we find that the lattice constants are underesti- 
mated by 1-2%. 



TABLE II. Relaxed structure of rhombohedral BaTi0 3 . 
Lattice constant a and atomic displacements A (relative to 
ideal cubic positions) in A; rhombohedral angle 8 in degrees. 





a 6 


A z (Ti) 


A 3,(0) 


A,(0) 


Theory 


4.00 89.85 


0.043 


-0.049 


-0.077 


Expt. a 


4.00 89.90 


0.052±12 


-0.044±8 


-0.072±8 



a Ref. 30. 
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TABLE III. Frequencies (in cm J ) of the zone-center opti- 
cal phonon modes in wurtzite ZnO. 



symmetry character 


Theory 


Expt. 31 


Expt 32 


Ai(TO) 


390 


380 


380 


Ai(LO) 


548 


574 


579 


Bi 


261 






Ei (TO) 


409 


407 


413 


Ei(LO) 


552 


583 


591 


E 2 


91 


101 


101 


E2 


440 


437 


444 



2. BaTiOz 

BaTi03 is a prototypical example of the class of per- 
ovskite ferroelectric materials. These materials normally 
have the paraelectric cubic perovskite structure at high 
temperature, but then undergo a ferroelectric instabil- 
ity as the temperature is reduced. BaTiC>3 actually goes 
through a series of three ferroelectric phase transitions as 
the symmetry is first tetragonal, then orthorhombic, and 
then rhombohcdral, with polarization respectively along 
[100], [110], and [111], with decreasing temperature. The 
ground-state rhombohedral structure (space group R3m, 
point group Cs v ) is fully determined by its lattice con- 
stant, rhombohedral angle, and the symmetry-allowed in- 
ternal atomic displacements along the [111] direction. We 
represent the rhombohedral phase in the hexagonal coor- 
dinate system, in which the z axis is along the previous 
[111] direction. Table II lists the structural parameters of 
our relaxed BaTiC>3, in which wc constrained the atomic 
volume to be equal to the experimental one as explained 
in Sec. III. The remaining structural parameters can be 
seen to be in good agreement with experiment. 

B. Displacement response tensors 

1. ZnO 

Wurtzite ZnO belongs to space group P63TOC {Cq v ). 
Standard group-theory analysis shows that the T-point 
phonon modes can be decomposed as 

r op t =Ax@ 2Bi © £1 © 2E 2 (40) 

in which the A\ and E\ modes are both Raman and 
IR active, while the nonpolar E 2 modes are Raman ac- 
tive only and B\ modes are silent. Shown in Table III 
are our computed phonon frequencies compared with two 
experimental results, showing good agreement with ex- 
periment. 

Because of the wurtzite symmetry of ZnO, the effec- 
tive charge tensor Z has only two independent elements, 
while the force-response internal-strain tensor A has only 



TABLE IV. Independent elements of the Born effective 
charge tensor Z (in units of e) and of the force-response in- 
ternal strain tensor (hartree/bohr) for wurtzite ZnO. 



Z xz {Zn) 


Z„(Zn) 






2.135 


2.163 






A x5 (Zn) 


A x6 (Zn) 


A z3 (Zn) 


A yl (0) 


-9.5 


-15.0 


18.7 


-16.7 



four independent elements. We present results for both 
tensors in Table IV. For this semiconductor material, it 
can be seen that the effective charge is very close to the 
nominal ionic charge. 

2. BaTiO A 

The low-temperature phase of BaTiOs has a rhombo- 
hcdral structure which belongs to the R3m space group. 
According to a group-theory analysis, the zone-center 
phonon frequencies can be decomposed as 

r op t - 3Ai © 4E © A 2 . (41) 

The Ai and E modes are both IR and Raman active, 
while the A 2 mode is silent. Table V gives the calcu- 
lated phonon frequencies at the T point. (The A 2 mode 
at 278 cm -1 and the E modes at 293 cm -1 are the ones 
derived from the silent F 2u modes of the undistorted cu- 
bic structure; because the distortion is weak, the LO-TO 
splitting of these E modes is negligible.) The results are 
very similar to those of the previous theoretical study 
of Ghosez. 33 While we are not aware of detailed experi- 
mental information about phonon frequencies at very low 
temperature, we note that measurements just below the 
orthorhombic to rhombohcdral phase transition temper- 
ature indicate phonon frequencies in three regions (100- 
300 cm -1 , 480-580 cm -1 , and 680-750 cm -1 ) in qualita- 
tive agreement with our zero-temperature calculations. 

We also calculated the atomic Born effective charges 
for this phase, but in view of the lower symmetry and 
larger number of independent elements, we have not 
listed them all here (our results are again very similar 



TABLE V. Phonon frequencies (in cm x ) at the T point 
for rhombohedral BaTiOs. 



Phonon mode 


Frequency 


Phonon mode 


Frequency 


Ai(TOl) 


169 


E(TOl) 


164 


Ai(LOl) 


179 


E(LOl) 


175 


Ai(T02) 


255 


E(T02) 


206 


Ai(L02) 


460 


E(L02) 


443 


Ai(T03) 


511 


E(T03) 


472 


Ai(L03) 


677 


E(L03) 


687 


A 2 


278 


E 


293 
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TABLE VI. Dielectric tensors of ZnO and BaTi0 3 (in units 
of e ). 







Present theory 


Experiment 




Index 




e e 


ZnO 


11 


5.76 10.31 11.09 


3.70 a 7.77 a 




33 


5.12 10.27 12.67 


3.78 a 8.91 a 


BaTi0 3 


11 


6.20 68.75 264.75 


6.19 b 




33 


5.79 37.44 49.51 


5.88 b 



a Ref. 34. 
b Ref. 35. 



to those of Ref. 33). The cation results are easily given 
as Z IX (Ba) = Z yy (Ba) = 2.78, Z zz (Ba) = 2.74, Z xx (Ti) 
= Z yy (Ti) = 6.64, and Z zz (Ti) = 5.83. The effective 
charge tensors are non-diagonal and non-symmetric for 
the oxygens; we mention only that the eigenvalues of the 
symmetric parts of these tensors cluster around —2 and 
—5, i.e., not much changed from their cubic-phase val- 
ues. Similarly, we have computed the full internal-strain 
tensor A for rhombohcdral BaTi03, but we have chosen 
not to present the details here because of the complicated 
form of this tensor involving a large number of indepen- 
dent elements. 



there is a soft ferroelectric mode present in the latter 
material. (Here, we use "soft" in the sense of a mode 
that has a small, but positive frequency; it is, of course, 
closely related to the imaginary-frequency unstable mode 
computed for the cubic structure, which condenses out to 
form the ferroelectric rhombohedral phase.) In the semi- 
conductor ZnO, on the other hand, no such soft mode is 
present. 

The last column presents our results for the free-stress 
dielectric tensors that are related to the fixed-strain 
tensors e via Eq. (20). The tensors e and are the 
same in higher-symmetry crystals, but in the presence of 
piezoelectric coupling, they are, in general, different. The 
free-stress tensors are always larger than the fixed-strain 
ones because the additional strain relaxation occurs so as 
to allow further polarization to develop in the direction of 
the applied field. We can see that the changes are modest 
for ZnO (on the order of 10-20%), which is not a soft- 
mode system. On the other hand, they are much more 
profound for the case of BaTiOs, where most notably an 
order-of-magnitude change occurs for en. This is related 
to the large value of the piezoelectric coupling factor /c 15 , 
as we will see later in Sec. II C 4. Essentially, it arises 
because the polarization vector is rather soft with respect 
to rotation away from the z axis, so that the electric 
susceptibility is large in the x-y plane. 



C. Dielectric tensors 

We now turn to a discussion of the computed dielec- 
tric tensors for wurtzite ZnO and rhombohedral BaTiOa, 
which are presented in Table VI. Because of the high 
point-group symmetry, the dielectric tensors have only 
two independent elements. Recall that the clampcd- 
ion tensor e, the fixed-strain relaxed-ion tensor e, and 
the free-stress relaxed-ion tensor are defined through 
Eqs. (14) and (17-20). While the results for the purely 
electronic dielectric tensors e are in good agreement with 
experiment for BaTiOa, we find that our LDA theory sig- 
nificantly overestimates the electronic dielectric response 
of ZnO. Hill and Waghmare, 22 also using an LDA pseu- 
dopotential approach, found €33 = 4.39, not as large as 
our 5.76, but still much larger than the experimental 
3.70. At least some of this overestimate is undoubtedly 
attributable to the LDA (and is connected with the un- 
derestimate of the gap in LDA), but the choice of pseu- 
dopotentials also seems to play a role. The computed 
lattice contributions en — en = 4.55 and €33 — 633 = 5.15 
arc in better agreement with the experimental values of 
4.07 and 5.13 respectively. 

While the clamped-ion tensors e are not so different 
for these two materials, the lattice contribution is clearly 
much bigger for the BaTi03 case. That is, while the lat- 
tice contribution (e — e) is about the same size as the 
purely electronic one (e) for ZnO, it is almost 10 times 
larger in BaTiOs. This difference clearly reflects the fact 



D. Elastic tensors 

We now consider the various elastic tensors. Recall 
that the clamped-ion elastic tensor C of Eq. (5) is just the 
second derivative of the unit-cell energy with respect to 
homogeneous strains, without allowing for internal struc- 
tural relaxations, while the physical elastic tensor C of 
Eq. (15) does include such relaxations. This tensor C 
(written more explicitly as C^') usually defined under 
conditions of fixed macroscopic electric field, but it is 
sometimes of interest to consider the elastic tensor 
of Eq. (23) defined instead under conditions of fixed elec- 
tric displacement field. These are identical for higher- 
symmetry (e.g., centrosymmetric) crystals, but that is 
not the case here. The compliance tensors S are defined 
as the inverses of the corresponding elastic tensors C. 

The results of our calculations of these tensors are dis- 
played in Table VII and VIII for ZnO and BaTiOs respec- 
tively. The lower point-group symmetry of BaTiOs is re- 
flected in the presence of an additional symmetry-allowed 
element C14. Actually, there are only five independent el- 
ements for ZnO, since 6*66 = (Cn — Ci2)/2 and Sqq — 
2(<Sii — Si 2) by symmetry. 17 Similarly, there are really 
only six independent elements for BaTiOs; in addition 
to the same relation, one also has Cqq = (Cn — C\^)/2, 
S 66 = 2(S n - S12), C 56 = C u , and S 56 = 2C U . Our 
results for the elastic constants of ZnO can be seen to be 
in good agreement with previous theory and with exper- 
iment (last columns of Table VII) . 
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TABLE VII. Clamped-ion (C) and relaxed-ion (C) elastic tensors at constant £, relaxed ion (C^ 13 - 1 ) elastic tensor at con- 
stant D (GPa), and corresponding compliance tensors (S, S, and S^) (TPa -1 ), for wurtzite ZnO. Previous theoretical and 
experimental results are also given for C for comparison. 



Index 


C 


C 




Present theory 
S 


S 


S iD) 


Theo. a 
C 


Expt. b 
C 


11 


305 


226 


231 


3.86 


7.79 


7.56 


246 


209 


12 


107 


139 


144 


-1.20 


-3.63 


-3.93 


127 


120 


13 


77 


123 


114 


-0.61 


-2.12 


-1.58 


105 


104 


33 


333 


242 


260 


3.29 


6.28 


5.23 


246 


211 


44 


62 


40 


43 


16.23 


24.69 


23.21 


56 


44 


66 


99 


44 


44 


10.12 


22.84 


22.73 







a Ref. 28. 
b Ref. 36. 



TABLE VIII. Clamped-ion (C) and relaxed-ion (C) elastic 
tensors at constant £, relaxed ion (C (D ^) elastic tensor at 
constant D (GPa), and corresponding compliance tensors (S, 
S, and S^ D ') (TPa- 1 ), for rhombohedral BaTiO a . 



Index 


C 


C 


C (D) 


S 


S 


S (D) 


11 


349 


277 


318 


3.32 


5.85 


3.65 


12 


106 


79 


93 


-0.82 


-2.94 


-0.95 


13 


96 


41 


81 


-0.72 


-0.45 


-0.68 


14 


8.4 


45 


19 


-0.31 


-8.17 


-0.89 


33 


334 


264 


323 


3.41 


3.93 


3.44 


44 


110 


48 


97 


9.12 


35.85 


10.63 


65 


8.3 


45 


19 


-0.63 


-16.33 


-1.78 


66 


121 


99 


113 


8.28 


17.58 


9.18 



We notice that the physical elastic Cjk are generally 
smaller than the frozen-ion ones Cjk (at least for diag- 
onal elements), since the additional internal relaxation 
allows some of the stress to be relieved. By the same 
token, diagonal S values are larger than S ones, reflect- 
ing the increased compliance allowed by the relaxation of 
the atomic coordinates. As for the dielectric constants, 
the differences are substantially smaller for ZnO than for 
BaTiC>3, as a result of the soft-mode contribution in the 
latter material. The constraint of fixed electric displace- 
ment field has the effect of suppressing some of this in- 
ternal relaxation (for essentially the same reason that 
longitudinal optical phonons are stiffer than transverse 
optical ones). This additional stiffness results in larger 
diagonal values than C values, and lower diagonal 

values than S values. However, the differences be- 
tween C^ and C tensors are generally smaller than the 
differences between C and C tensors, especially for ZnO. 



E. Piezoelectric tensors 

The bare (or "frozen-ion") piezoelectric tensor e a j is 
just given by the mixed second derivatives of unit cell 
energy with respect to electric field £ a and strain r]j , de- 
forming internal atomic coordinates in strict proportion 
to the homogeneous strain. The full piezoelectric tensor 
e a j also takes into account the contributions from the lat- 
tice, as described in Eq. 16. The total number of indepen- 
dent piezoelectric tensor members is determined by the 
point group of material. Rhombohedral BaTiOa (point 
group C3,,) has a lower symmetry than that of wurtzite 
ZnO (Cq v ), so we may expect more independent elements 
in the former. Indeed, a symmetry analysis 17 shows that 
ZnO has only three independent tensor elements, namely 
e3i and 633 describing polarization along the c axis in- 
duced by uniaxial c-axis or biaxial a6-plane strains, while 
e24 describes the polarization induced by shear strains. 
For BaTiOs, the symmetry is slightly lower, and as a re- 
sult there is a fourth independent tensor element in this 
case. 

In Table IX and X, we present our results for piezo- 
electric tensors for these two materials. We also also give 
the results for the d a j tensor as defined in Eq. 30. Our 
results for the e a j matrix for ZnO are consistent with 
the previous theory. 22 ' 37 ' 38 (Table X shows five tensor el- 
ements for BaTiOs, not four, but in fact e2i = e±e and 
2g?2i = c£i6 by symmetry. 17 ) 

Recall that the frozen-ion and relaxed-ion piezoelec- 
tric tensors are defined by Eqs. (8) and (16), in which 
the relaxed-ion piezoelectric tensor incorporates contri- 
butions from lattice relaxation. For the same reason as 
discussed previously for the case of the dielectric and elas- 
tic tensors, the difference between the above two tensors 
(e.g. e and e) is much bigger for BaTiOs than for ZnO, 
as expected from the presence of the soft mode in the 
perovskite material. Also, note that the electronic and 
lattice contributions have opposite signs, with the lattice 
contribution being the larger of the two, as is common 
for other dielectric materials. 
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TABLE IX. Clamped-ion e (C/m 2 ), relaxed-ion e (C/m 2 ), and relaxed-ion d (pC/N) piezoelectric tensors for wurtzite ZnO. 
Previous theoretical and experimental results are also given for e and d for comparison. 



Index 


e 


Present theory 

e 


d 


Theo. a 

e 


d 


Expt. b 

e 


d 


31 


0.37 


-0.67 


-5.5 


-0.55 


-3.7 


-0.62 


-5.1 


33 


-0.75 


1.28 


10.9 


1.19 


8.0 


0.96 


12.3 


15 


0.39 


-0.53 


-13.1 


-0.46 


-8.2 


-0.37 


-8.3 



a Ref. 28. 
b Ref. 36. 



TABLE X. Clamped-ion e (C/m 2 ), relaxed-ion e (C/m 2 ), 
and relaxed-ion d (pC/N) piezoelectric tensors for rhombohe- 
dral BaTi0 3 . 



Index e e d 

21 -0.23 2lH 70T 

31 0.05 -3.03 -6.8 

33 -0.19 -4.44 -14.7 

15 0.18 -5.45 -243.2 

16 -0.23 2.91 140.2 



In view of this partial cancellation of terms of opposite 
sign, accurate calculations of e and d coefficients are par- 
ticularly delicate. Wc find that our results for the 31 and 
33 elements of the e and d coefficients of ZnO are in rea- 
sonably good agreement with experiment (slightly better 
than previous Hartree-Fock calculations 28 ), whereas we 
somewhat overestimate the shear coefficients eis and d\§ 
(slightly more so than in the Hartree-Fock theory 28 ). 

F. Electromechanical coupling constants 

We compute and present in Table XI the piezoelectric 
coupling factors /C33, £31, and fcis defined in Eq. (34) for 
both ZnO and BaTiOs. We also calculate the singular 
values fc„ of the /C matrix of Eq. (35) . Because of the axial 
symmetry, these are arranged into a pair of degenerate 
values ki = corresponding to in-plane electric fields, 
and a fc 3 corresponding to axial fields. (In fact, due to 
the symmetry, fc„ cab be calculated in practice as just 
[f3il> (d ■ ■ d T ) uv ] 1 /' 2 .) These k v values are also given 
in Table XI. 

Roughly speaking, the couplings given in the first 
three lines of Table XI are associated with symmetry- 
preserving "polarization stretching" degrees of freedom, 
while those in the last two lines correspond to "polariza- 
tion rotation" modes. Note that ki 5 = k\ for ZnO but 
not for BaTiOs; this is a feature of symmetry, arising be- 
cause an electric field Z\ couples uniquely to 775 in ZnO, 
but ajso to r/6 in BaTi03. Also, we can see that k\ > fcis 
and ki > max(fc33,fc3i) in both materials, since k de- 
scribes the optimal coupling between electric and elastic 



TABLE XL Dimensionless piezoelectric coupling factors. 
The first three correspond to £ -fields longitudinal to the crys- 
tal axis; the last two are transverse. Coupling constants fci 
and k 3 are obtained from singular-value analysis of the cou- 
pling tensor K. (see text). 



ZnO BaTiOs 



k 33 0.41 0.35 

fc 3 i 0.19 0.13 

k 3 0.44 0.49 

fcis 0.27 0.84 

fei 0.27 0.86 



channels. 

Comparing the two materials, we see that the coupling 
factors are rather comparable in the polarization stretch- 
ing channel; evidently, the soft mode does not play such a 
profound role there. In contrast, the coupling factor fci 5 
is very large in BaTiOs; in fact, it is not far from unity, 
the maximum value consistent with stability. Indeed, this 
is precisely because the crystal is not far from being un- 
stable with respect to a rotation of the polarization away 
from the rhombohcdral axis - precisely the type of dis- 
tortion that would carry it to the orthorhombic phase, 
from which it evolved as the temperature was reduced 
below the orthorhombic-rhombohcdral phase transition 
temperature that occurs experimentally at ~ 180 K. The 
large fcis is also strongly connected to the large difference 
between e = and in Tab. VI as already discussed 
at the end of Sec. IV C. 

The calculated value fci5 = 0.84 for BaTiOs is quite 
respectable; it is in the range of the values of fcis = 0.25- 
0.80 39 for the PMN-PT and PZN-PT single-crystal piezo- 
electrics on the rhombohcdral side of the morphotropic 
phase boundary. Unfortunately, the fact that this large 
coupling occurs only at very low temperature probably 
makes it useless for practical applications. On the other 
hand, the present work suggests that if a material like 
BaTiOs could somehow be stabilized in the rhombo- 
hcdral phase at room temperature, it might have very 
promising piezoelectric properties. 
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V. SUMMARY 

In summary, we have developed a method that system- 
atically treats the effects of perturbations associated with 
atomic displacements, electric fields and strains in insu- 
lating crystals, so that physical quantities expressible as 
second derivatives of the total energy can be computed 
efficiently In the first step, six elementary tensors are 
computed once and for all using the methods of DFPT: 
the force-constant matrix, the Born charge tensor, the 
internal-strain tensor, and the frozen-ion dielectric, elas- 
tic, and piezoelectric tensors. The internal-displacement 
degrees of freedom are then eliminated to give physical 
low-frequency dielectric, elastic, and piezoelectric ten- 
sors, defined under boundary conditions of controlled 
electric field £ and strain rj. We have also shown how 
these can then be manipulated to obtain various related 
tensor properties of interest such as the free-stress dielec- 
tric tensor, the fixed-D elastic and compliance tensors, 
and various piezoelectric tensors and electromechanical 
coupling factors. Such a systematic approach is espe- 
cially important in polar crystals, in which the atomic- 
displacement, electric-field, and strain degrees of freedom 
are strongly coupled in complex ways. 

We have applied our approach to compute these ten- 
sor properties for two paradigmatic crystals, ZnO and 
BaTi03, at zero temperature. These materials differ 
most significantly in that there is a ferroelectric soft mode 
that has condensed, but still remains rather low in fre- 
quency, in the latter material. The calculations are sub- 
ject to several approximations, most notably the LDA 
itself (and its associated lattice-constant error, which 
has been removed by hand for the case of BaTiC>3 - see 
Sec. Ill), but also the frozen-core approximation (as im- 
plemented through the use of pseudopotentials) and the 
neglect of zero-point fluctuations. Nevertheless, we vali- 
date the approach by finding reasonably good agreement 
between theory and experiment for most quantities in 
the case of ZnO, despite the fact that the experiments are 
room-temperature ones. The largest discrepancies are for 
the purely electronic dielectric tensor elements en and 
£33, the shear piezoelectric coupling ei5, and to a lesser 
extent, derived quantities that depend on these elemen- 
tary ones. In the case of BaTiC>3, where low-temperature 
experiments on single-crystal, single-domain samples un- 
der well-defined boundary conditions are not available, 
our calculations provide useful predictions of the material 
constants. In particular, we find an encouraging value of 
0.84 for the kis electromechanical coupling constant, and 
argue that this is associated with the proximity of the or- 
thorhombic phase. 

We wish to emphasize that the usefulness of the general 
approach advocated here transcends the particular imple- 
mentation of it (here DFT/LDA, pseudopotentials, etc.). 
For example, similar calculations might by carried out 
with Hartree-Fock methods using localized orbitals 28 or, 
eventually, using "LDA+U", dynamical mean-field the- 



ory, or quantum Monte Carlo methods. In this case, the 
six elementary tensors of Eqs. (3-8) will first need to be 
calculated using methods appropriate to the particular 
type of electronic-structure method used. However, the 
subsequent manipulations described in Sees. II B and II C 
can then be carried through in identically the same way 
as done here. 

Finally, we note that the approach described here can 
be extended to include other types of perturbations, such 
as alchemical ones, and to the treatment of higher-order 
responses (e.g., anharmonic elastic constants, nonlinear 
dielectric responses, and electrostriction effects), provid- 
ing possible directions for future developments of the 
method. 
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APPENDIX A: SIMULTANEOUS TREATMENT 
OF STRAINS AND ELECTRIC FIELDS 

The formulation of an energy functional, and the defi- 
nition of response functions in terms of its second deriva- 
tives, is somewhat subtle in the case that electric fields 
and strains are simultaneously present. The purpose of 
this Appendix it to give a careful treatment of the theory 
in this case. Except where noted, the notation here fol- 
lows that of Sec. II B in that the internal displacements 
assumed to have been integrated out (i.e., internal 
displacements u m do not appear explicitly). 

We begin introducing the deformation tensor rj a f3 in 
the Cartesian frame via 



dr 



(Al) 



(implied sum notation) where dr a is the deformation 
of the medium from its undeformed position r a . We 
consider deformations taking the form of homogeneous 
strains and rigid rotations, so that the antisymmetry part 
of the deformation tensor rj a p describes the rotational 
part, while its symmetric part is just the strain tensor 



1 ,~ 

The improper piezoelectric tensor is defined as 

dP a 



o impr 

'a/37 



(A2) 



(A3) 



The name "improper" reflects that fact that e lmpr con- 
tains contributions that are spurious in a certain sense. 15 
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For example, if we consider a pure rotation of a spon- 
taneously polarized crystal about an axis that does not 
coincide with P, then dP does not vanish, and conse- 
quently e 1 ™^* has a component that is antisymmetric un- 
der interchange of and 7. Similarly, if we consider a 
uniaxial or biaxial compression in the x — y plane of a 
ferroelectric modeled as a lattice of discrete rigid dipoles 
oriented in the z direction, the polarization will change 
even if the dipoles do not, because the polarization is de- 
fined as the dipole moment per unit volume. This, too, is 
an essentially spurious effect. By contrast, the "proper" 
piezoelectric tensor will be defined so as to vanish in ei- 
ther of these situations. 

In the presence of a strain, it is convenient to introduce 
reduced or rescaled polarizations P' a and electric fields £' a 
via 



and 



£' a = rjfj a £p 



(A4) 



(A5) 



where P' a and £' a coincide with P a and £ a in the absence 
of strains or rotations. We then take our fundamental 
energy functional to be 



"0 



Ei%-Q£'.P' 



(A6) 



where is the usual zero-field Kohn-Sham energy per 
cell 14 of the occupied Bloch functions introduced earlier 
in Sec. II A. 

Note that P • £ = P' • £ ', so that Eq. (A6) is closely re- 
lated to Eq. (1). However, it is important to understand 
that £' and rj are the "natural variables" of the energy 
functional (A6), so that subsequent partial derivatives 
are defined in terms of this pair of variables. For exam- 
ple, the proper piezoelectric tensor is now given by 



d 2 E 



d£' a dr/p-, 



or 



1 d{np' a ) 

Q dripj 



(A7) 



(A8) 



We also emphasize that £' is, in many ways, a more 
natural variable than £ from the experimental point of 
view. For example, if one controlls the voltage V across 
a film of M atomic layers between conducting capacitor 
plates and observes the resulting strain, one is actually 
controlling £' = eV/Mc n = £c/c n , not £, where cq and 
c are the zero-field and finite-field lattice constants, re- 
spectively, in the normal direction. 

From Eqs. (A3) and (A8) it follows that the improper 
and proper piezoelectric tensors are related by 15 



e a/ 3 7 = e^ r - 5p 1 P a + S a fj Pj 



It is then easy to show that the proper tensor e a py is 
symmetric under interchange of indices (3 and 7, so that 
the Voigt notation can be restored. This is to be expected 
because the reduced quantity P' is invariant under a rigid 
rotation of the crystal, a fact that follows trivially from 
Eq. (A4). That is, P', expressed as a function of the six 
symmetrized strain variables (Eq. (A2)) and the three 
rotational variables, is actually independent of the rota- 
tional ones. Indeed, a rigid rotation of the entire sys- 
tem, material plus external field £, leaves both £' and 
P' individually unchanged. It is then natural to discard 
the rotational variables and recast the symmetric strain 
variables in Voigt notation. We then regard the energy 
functional of Eq. (A6) to be a functional E(£', 77) of the 
fundamental variables of rescaled £' a and Voigt rjj, and 
the proper piezoelectric tensor may be written as 



e a j — 



d 2 E 1 d(np^) 



d£' a drjj Q drjj 



(A10) 



Restoring the explicit dependence on internal displace- 
ments u m , Eq. (A6) becomes 



E{u,£',ri) 



1 



E^-n£>.p> 



(All) 



where E(£' ,rj) of Eq. (A6) corresponds to the minimum 
of (All) over all possible displacements u m . While this 
equation is numerically equal to Eq. (1), it is critical to 
recall that it is written in terms of different arguments. 

Strictly speaking, this notation should have been in- 
troduced at the very beginning of Sec. II A, and every 
equation throughout the paper, starting with Eq. (1), 
should have £ replaced by £' and P by P'. For example, 
Eqs. (8) and Eq. (26) should be replaced by Eqs. (A12) 
and 



d 2 E 

d£' a d Vj 



(A12) 



respectively, and similarly for all other equations. How- 
ever, for the purposes of clarity of presentation, it was 
decided to avoid use of this clumsy notation in the main 
part of the paper. 

Finally, we note that the reduced quantities P' and £' 
are also rather natural physical variables from the point 
of view of computational implementation. Indeed, these 
two quantities can further be expressed in terms of fully 
reduced quantities p^ and e M via 



(0) 



and 



P' — — n a' 



e£'-a<°> 



(A13) 



(A14) 



so that P = (e/fi)p M a M and e M = e£ ■ a M . In these equa- 
(A9) tions, a^ is the ^i'th undeformed primitive real-space 
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lattice vector, a M is the corresponding deformed lattice 
vector, and £1 — ai • a 2 x a 3 is the deformed cell volume. 
Note that P • £ = P' • £' = ft" 1 p M e M and 



fin 



(A15) 



The fully reduced polarization p M has a simple interpre- 
tation in terms of the fractional positions of the charges 
in the unit cell; for example, the contribution to it com- 
ing from filled band is just — l/2n times the Berry phase 
of that band, as can be seen by comparing with Eq. (10) 
of Ref. 15. Similarly, is just e times the electrostatic 
potential drop across the unit cell in direction a M . 

The computational implementation of DFPT is done 
quite naturally in terms of these reduced quantities, 7 ' 10,19 
and as a result, DFPT automatically yields the proper 
piezoelectric tensor. 10 This can be a source of confu- 
sion when comparing the DFPT results with those of 
finite-difference calculations. In the latter approach, 
the polarization is obtained directly from ground-state 
DFT calculations, 6 and piezoelectric tensor elements are 
obtained by numerical differentiation using sufficiently 
small strains about the reference structure. This proce- 
dure yields the improper tensor, however, and Eq. (A9) 
must be applied to compare such results with the DFPT 



ones 
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